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Abstract 

Reduced  order  modeling  (ROM)  seeks  to  systematically  extract  the  important  dynamics  from  a 
high  fidelity,  high  dimensional  dynamical  system  such  that  outputs  of  interest  generated  by  the 
high  fidelity  simulation  and  the  ROM  over  a  range  of  input  parameters  are  within  user-specified 
bounds.  ROM  is  important  in  real-time  applications,  design,  control,  and  probabilistic  analyses, 
where  outputs  of  interest  need  to  be  evaluated  quickly  for  various  input  parameters. 

This  collaborative  research  has  developed  rigorous  mathematical  and  computational  frame¬ 
works  for  ROM  generation  and  the  use  of  ROMs  in  real-time,  design,  control,  and  probabilistic 
applications  of  relevance  to  the  Air  Force.  This  research  has  provided  theoretical  analyses  and  nu¬ 
merical  studies  for  several  new/extensions  of  existing  ROM  approaches,  such  as  a  goal-oriented, 
model-constrained  approach,  balanced  truncation  model  reduction  (BTMR)  of  descriptor  systems, 
and  the  integration  of  domain  decomposition  and  BTMR  for  systems  with  localized  nonlinearities. 
Additionally,  several  important  questions  related  to  the  design,  analysis,  efficient  computation,  and 
application  of  ROM  were  studied.  The  use  of  ROMs  was  demonstrated  on  example  applications, 
including  optimal  flow  control  of  linearized  Navier-Stokes  equations,  linearized  flow  control  of  a 
supersonic  diffuser,  subsonic  compressor  blade  row  unsteady  aerodynamics  and  geometric  mistun- 
ing,  a  thermal  design  problem,  and  nonlinear  combustor  model  dynamics. 
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1  Introduction 


With  the  advent  of  efficient  algorithms  and  high  performance  computers,  the  use  of  high-fidelity 
numerical  simulations  of  fluid  systems  is  moving  from  single-run  calculations  that  explore  flow 
physics  towards  multi-run  engineering  optimization  calculations  and  probabilistic  analyses.  These 
studies  require  the  repeated  numerical  simulation  of  a  complex  system  to  determine  outputs  of 
interest  for  varying  design/control/input  parameters. 

Reduced  order  modeling  seeks  to  systematically  extract  the  important  dynamics  from  a  high 
fidelity  dynamical  system,  involving  a  huge  number  of  degrees  of  freedom,  in  such  a  way  that  the 
error  between  outputs  of  interest  generated  by  the  costly  high  fidelity  simulation  and  the  inexpen¬ 
sive  low  dimensional  reduced  order  model  (ROM)  over  a  range  of  input  parameters  is  below  a 
user-specified  tolerance. 

Although  there  has  been  gratifying  progress  in  design  and  analysis  of  ROM  methods  (see, 
e.g.,  [3])  and  the  application  of  ROM  to  aerospace  systems  (see  e.g.,  [16]),  the  design  of  ROM 
approaches,  derivation  of  error  bounds,  and  the  efficient  computational  implementation  of  ROM 
approaches  for  systems  with  relevance  to  Air  Force  applications  pose  many  challenges.  This  col¬ 
laborative  research  has  successfully  tackled  several  important  theoretical  and  algorithmic  issues 
related  to  the  design,  analysis,  and  efficient  computation  of  ROM  for  problem  classes  which  are 
important  in  Air  Force  applications.  Furthermore,  the  use  of  ROMs  was  demonstrated  on  several 
example  applications  with  relevance  to  the  Air  Force,  including  optimal  flow  control  of  linearized 
Navier-Stokes  equations,  linearized  flow  control  of  a  supersonic  diffuser,  subsonic  compressor 
blade  row  unsteady  aerodynamics  and  geometric  mistuning,  a  thermal  design  problem,  and  nonlin¬ 
ear  combustor  model  dynamics.  The  techniques  developed  in  this  research  have  also  been  applied 
to  reduced  order  modeling  of  nonlinear  systems  arising  in  the  computational  simulation  of  large 
neuron  systems. 


2  Accomplishments 

This  project  has  considered  projection  based  model  reduction  methods.  Abstractly  the  original 
system  (typically  a  discretized  system  of  partial  differential  equations)  can  be  written  as 

M=/(y(*W)  *€(o,r),  (ia) 

z(t)=g(y(t),u,t)  te(0,T),  (lb) 

where  y(t)  G  RA’  is  the  system  state  (y(t)  are,  for  example,  velocities  and  pressure  at  time  t),  u  are 
the  possibly  time  dependent  inputs  (design  variables,  controls,  uncertain  system  parameters),  and 
z(t )  G  R*  are  the  outputs  of  interest.  We  seek  a  reduced  order  model  of  he  form 

%)  =  <rtTf{Vy{t),u,t)  t  G  (0,r),  (2a) 

m  =g(Vy{t)(t),u,t)  re(0,r),  (2b) 

where  y(t)  G  R”  is  the  reduced  system  state,  n  <^N.  As  before,  u  are  the  inputs  and  z(t)  G  R*  are 
the  outputs  of  interest,  which  are  now  generated  via  a  reduced  state  equation  (2a).  The  matrices 
G  R^x”  satisfy  tWT(U  =  /. 
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This  research  has  addressed  several  important  questions  related  to  the  computation  of  the  pro¬ 
jection  matrices  V,  e  RNxn,  estimates  for  the  error  ||z  —  z||  between  outputs  of  the  full  and  the 
reduced  order  model,  and  fast  computation  of  reduced  order  models  for  nonlinear  models.  Specif¬ 
ically,  this  research  investigated  a  goal-oriented,  model-constrained  approach,  balanced  truncation 
model  reduction  (BTMR)  of  descriptor  systems,  and  the  integration  of  domain  decomposition  and 
BTMR  for  systems  with  localized  nonlinearities.  Several  important  questions  related  to  the  design, 
analysis,  efficient  computation,  and  application  of  ROM  in  the  context  of  applications  of  relevance 
to  Air  Force  were  studied. 

Furthermore,  the  use  of  ROMs  was  demonstrated  on  several  test  problems  which  are  relevant 
for  Air  Force  applications,  including  optimal  flow  control  of  linearized  Navier-Stokes  equations, 
linearized  flow  control  of  a  supersonic  diffuser,  subsonic  compressor  blade  row  unsteady  aero¬ 
dynamics  and  geometric  mistuning,  a  thermal  design  problem,  and  nonlinear  combustor  model 
dynamics.  The  techniques  developed  in  this  research  have  also  been  applied  to  reduced  order 
modeling  of  nonlinear  systems  arising  in  the  computational  simulation  of  large  neuron  systems. 

The  findings  of  this  research  are  described  below.  Additional  details  and  references  may  be 
found  in  the  papers  [2,  5,  6,  7,  8,  12,  13,  14,  15,  18,  19]  and  theses  [4,  9,  10,  17]. 

2.1  Balanced  Truncation  Model  Reduction  for  the  Linearized  Navier-Stokes 
Equations 

Balanced  Truncation  Model  Reduction  (BTMR)  is  a  particular  model  reduction  method  that  pre¬ 
serves  asymptotic  stability  and  also  provides  an  error  bound  on  the  discrepancy  between  the  outputs 
of  the  full  and  reduced  order  system.  Unfortunately,  the  vast  majority  of  papers  on  BTMR  apply 
this  technique  to  dynamical  systems  of  ordinary  differential  equations,  but  are  not  directly  appli¬ 
cable  to  dynamical  systems  governed  by  incompressible  flows,  which  lead  to  large-scale  system  of 
differential  algebraic  equations. 
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Figure  1 :  Left  plot:  Time  response  for  the  full  order  model  (circles)  and  for  the  reduced  order  model 
(solid  line).  Right  plot:  Velocities  generated  with  the  22K  dof  full  order  model  (left  column)  and 
with  the  15  dof  reduced  order  model  (right  column)  at  t  =  1,3,4, 5,6  (top  to  bottom). 


The  work  [12, 17]  describes  an  extend  on  of  BTMR  to  a  class  of  differential  algebraic  dynamical 
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systems  including,  e.g.,  linearized  Navier-Stokes  equations,  linearized  around  steady  flow.  The 
procedure  in  [12, 17]  can  serve  as  a  template  for  handling  numerous  other  similar  problems  arising 
from  CFD.  The  major  advantages  of  this  approach  are  that  it  operates  on  the  original  equations 
and  not  the  projected  ones,  that  it  considerably  reduces  storage  requirements,  and  that  it  utilizes 
more  straightforward  linear  algebra  techniques  based  upon  saddle  point  solvers  that  are  most  likely 
already  available  to  a  user.  Finally,  the  approach  in  [12,  17]  produces  reduced  order  models  with 
guaranteed  error  bounds. 

To  illustrate  the  last  point,  the  linearized  Navier-Stokes  equations  in  a  backward  facing  step 
domain  is  considered.  The  inputs  are  suction/blowing  controls  applied  on  the  inflow  and  backward 
facing  step  boundary,  the  output  of  interest  is  the  integrated  vorticity  in  the  region  behind  the 
backward  facing  step.  Figure  1  shows  an  excellent  agreement  between  the  time  response  of  the 
full  order  model  and  of  the  reduced  order  model  for  an  arbitrary  input.  This  excellent  agreement 
is  a  consequence  of  the  theoretical  error  bound  between  the  original  and  reduced  order  model. 
The  right  plot  Figure  1  shows  the  velocity  profiles  generated  by  the  original  and  the  reduced  order 
system.  There  are  slight  differences,  but  these  are  to  be  expected  since  the  balanced  reduction  was 
designed  generate  a  reduced  order  model  that  faithfully  reproduces  the  integral  of  vorticity. 

2.2  Numerical  Solution  of  Large  Scale  Lyapunov  Equations  Using  Inexact 
Linear  System  Solves 

A  computationally  expensive  part  of  BTMR  is  the  computation  of  approximate  controllability  and 
observability  Gramians.  If  one  computes  them  as  approximate  solutions  of  the  controllability  and 
observability  Lyapunov  equations  via  iterative  methods,  such  as  the  multishift  ADI  Algorithm,  then 
in  each  iteration  one  has  to  solve  large  scale  linear  systems  resulting  from  PDE  discretizations.  For 
very  large  scale  problems,  especially  3D  problems,  these  systems  have  to  be  solved  iteratively.  The 
thesis  [17]  gives  a  detailed  analysis  of  the  impact  of  these  iterative,  hence  inexact  linear  system 
solves  on  the  convergence  of  the  ADI  Algorithm  as  well  as  on  the  quality  of  the  computed  solution. 
Inexact  linear  system  solves  have  little  impact  on  the  convergence  behavior  as  well  as  on  the  quality 
of  the  computed  solution,  provided  the  stopping  criteria  for  the  linear  system  solves  are  sufficiently 
fine  relative  to  the  ADI  convergence  rate  and  stopping  tolerance. 


2.3  Coupling  with  Local  Nonlinearities 

BTMR  is  one  of  the  very  few  model  reduction  techniques  for  which  error  bounds  between  the 
input-output  map  of  the  full  and  reduced  order  model  are  available.  Unfortunately,  these  are  only 
available  for  linear  time  invariant  systems.  In  this  research  we  have  extended  BTMR  to  systems 
with  spatially  localized  nonlinearities  [17,  18].  Domain  decomposition  techniques  are  used  to 
divide  the  problem  into  Unear  subproblems  and  small  nonUnear  subproblems.  It  also  identifies 
interface  conditions  between  the  Unear  and  the  nonUnear  subdomains.  BTMR  is  appUed  to  the 
Unear  subproblems  with  inputs  and  outputs  determined  by  the  original  in-  and  outputs  as  weU  as 
the  interface  conditions  between  the  subproblems.  Figure  2  shows  the  result  of  our  approach  for 
a  PDE  that  couples  the  nonUnear  Burgers  equation  on  [—1, 1]  to  the  heat  equation  on  [-10,-1], 
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[1,10].  The  linear  parts  are  reduced  and  joined  with  the  nonlinear  Burgers  equation.  The  plots 
show  an  excellent  agreement  between  full  and  reduced  order  models. 


Figure  2:  Left:  Solution  of  the  reduced  order  discretized  PDE.  Center:  Error  between  the  solution 
of  the  discretized  PDE  and  the  reduced  order  system.  Right:  Outputs  1 , 2, 3  of  the  full  order  system 
are  given  by  *,  o  and  □.  Outputs  1,  2,  3  of  the  reduced  order  system  are  given  by  dotted,  dashed 
and  solid  lines. 

The  thesis  contains  [17]  contains  additional  results  for  other  test  problems. 

2.4  Model  Reduction  for  Large-Scale  Systems  with  High-Dimensional  Para¬ 
metric  Input  Space 

The  computation  of  projection  matrices  V,  TV  G  RNxn,  such  that  the  error  ||z- z]|  between  out¬ 
puts  of  the  full  and  the  reduced  order  model  is  minimized  over  a  range  of  system  inputs  u  can 
be  viewed  as  an  optimization  problem.  The  approximate  solution  of  this  optimization  problem 
is  very  challenging,  especially  for  parametric  input  spaces  of  high  dimension.  We  have  proposed 
a  model-constrained  adaptive  sampling  methodology  for  reduction  of  large-scale  systems  with 
high-dimensional  parametric  input  spaces.  Our  model  reduction  method  uses  a  reduced  basis  ap¬ 
proach,  which  requires  the  computation  of  high-fidelity  solutions  at  a  number  of  sample  points 
throughout  the  parametric  input  space.  The  papers  [4,  8]  pose  the  task  of  determining  appro¬ 
priate  sample  points  as  a  partial  differential  equation  (PDE)  -  constrained  optimization  problem, 
which  is  implemented  using  an  efficient  adaptive  algorithm  that  scales  well  to  systems  with  a  large 
number  of  parameters.  For  a  heat  transfer  optimal  design  application  with  18,000  states,  the  ap¬ 
proach  is  demonstrated  for  parametric  input  spaces  up  to  dimension  21  (see  [8]).  Sampling  such 
high-dimensional  input  spaces  with  statistically-based  sampling  methods  is,  in  general,  a  computa¬ 
tionally  prohibitive  proposition;  the  model-constrained  sampling  developed  as  pan  of  this  research 
yields  reduced  models  with  errors  three  to  four  orders  of  magnitude  lower  than  those  obtained 
using  standard  methods,  such  as  Latin  hypercube  and  log-random  sampling. 
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CFD 

Reduced 

Model  size 

103,008 

201 

Offline  cost 

- 

2.8  hours 

Online 

501.1  hours 

0.2 1  hours 

Blade  1  mean 

-1.8572 

-1.8573 

Blade  1  variance 

2.687e-4 

2.682e-4 

Blade  2  mean 

-1.8581 

-1.8580 

Blade  2  variance 

2.797e-4 

2.799e-4 

Table  1:  Linearized  CFD  model  and  reduced-order  model  MCS  results.  Work  per  cycle  is  predicted 
for  blade  plunging  motion  at  an  interblade  phase  angle  of  180°  for  10,000  randomly  selected  blade 
geometries. 

The  method  has  been  applied  to  derive  efficient  reduced-order  models  for  probabilistic  analysis 
of  the  effects  of  blade  geometry  variation  for  a  two-dimensional  model  problem  governed  by  the 
Euler  equations  [5,  6].  The  CFD  model  uses  a  discontinuous  Galerkin  formulation,  and  has  51,504 
states  per  blade  passage.  Figure  3  shows  the  results  of  a  Monte  Carlo  simulation  (MCS)  to  analyze 
the  impact  of  blade  geometry  variabilities  on  the  work  per  cycle  (WPC),  which  is  defined  as  the 
integral  of  blade  motion  times  the  lift  force  over  one  unsteady  cycle.  The  case  shown  considers 
two  blades  moving  in  plunging  motion  with  an  interblade  phase  angle  of  180°.  The  same  10,000 
geometries  were  analyzed  using  the  CFD  model  and  a  reduced  model  of  dimension  201.  Figure  3 
shows  the  resulting  probability  density  functions  of  WPC  for  the  first  blade.  Table  1  shows  that  the 
CPU  time  required  to  compute  the  reduced  model  MCS  is  a  factor  of  more  than  2000  times  smaller 
than  that  required  for  the  CFD  MCS.  Table  1  also  compares  the  statistics  of  the  two  distributions.  It 
can  be  seen  that  the  reduced-order  model  predicts  the  mean,  variance  and  shape  of  the  distribution 
of  WPC  accurately.  To  further  verify  the  quality  of  the  reduced  model,  the  Kolmogorov-Smimov 
method  is  applied  to  test  whether  the  reduced  WPC  results  and  the  full  WPC  results  are  drawn 
from  a  same  distribution.  The  results  show  that  the  hypothesis  that  the  distribution  is  the  same  at  a 
5%  significance  level  cannot  be  rejected. 
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wpc 

(a)  Full  WPC  for  Blade  1 


wpc 

(b)  Reduced  WPC  for  Blade  1 


Figure  3:  Comparison  of  linearized  CFD  and  reduced-order  model  predictions  of  WPC  for  Blade 
1.  MCS  results  are  shown  for  10,000  blade  geometries.  The  same  geometries  were  analyzed  in 
each  case. 

2.5  Efficient  Evaluation  of  Projection  Based  Reduced  Order  Models  for  Non¬ 
linear  Systems 

Methods  such  as  the  proper  orthogonal  decomposition  (POD)  and  the  optimization-based  approach 
[5,  6].  can  be  used  to  compute  a  reduced  basis  V  =  W  for  nonlinear  systems;  however,  naive 
projection  of  the  governing  equations  onto  the  reduced  basis  can  lead  to  a  reduced  model  that  has  a 
small  number  of  states  but  is  computationally  expensive  to  solve,  since  evaluation  of  the  nonlinear 
term 

y(0. »-»  *->  /( Vy(t) >  0,  ■->  yr/( 0 

eRn  eRN  eRN  eR" 

will  in  general  require  computations  on  the  large  scale  (R,v). 

This  research  has  combined  the  model-constrained  sampling  algorithm  described  in  Section  2.4 
with  the  Empirical  Interpolation  Method  of  [1,  1 1]  for  efficient  representation  of  nonlinearities.  In 
this  approach,  the  nonlinearities  are  represented  using  a  coefficient-function  approximation  that 
enables  the  development  of  an  efficient  offline-online  computational  procedure  where  the  online 
computational  cost  is  independent  of  the  size  of  the  original  large-scale  model. 
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Figure  4:  Comparison  of  (a)  finite  element  solution  and  (b)  reconstructed  solution  using  the  effi¬ 
cient  reduced  basis  approximation  for  nonlinear  convection-diffusion-reaction  problem. 

The  nonlinear  model  reduction  methodology  was  applied  to  a  highly  nonlinear  combustion 
problem  governed  by  a  convection-diffusion-reaction  PDE  [10].  The  reduced  basis  approximation 
developed  for  this  problem  is  up  to  50,000  times  faster  to  solve  than  the  original  high-fidelity 
finite  element  model  with  an  average  relative  error  in  prediction  of  outputs  of  interest  of  2.5  •  10-6 
over  the  input  parameter  space.  Figure  4  shows  a  sample  comparison  between  the  solution  field 
computed  using  the  efficient  reduced  basis  approximation  and  the  truth  finite  element  solution  for 
a  random  point  in  the  parameter  space.  The  reconstructed  field  is  virtually  indistinguishable  from 
the  original  finite  element  solution. 

Additionally,  an  alternative  view  of  the  Empirical  Interpolation  Method  was  developed  in  [9].  It 
is  applicable  to  fully  discretized  systems  and  uses  matrix-vector  representations  of  the  discretized 
system,  rather  than  the  variational  form  of  the  PDE  and  Galerkin  discretization.  The  presentation 
in  [9]  also  provides  alternative  error  bounds  for  EIM  based  on  matrix  analysis. 

Additionally,  the  nonlinear  model  reduction  based  on  POD  and  EIM  was  applied  to  a  nonlinear 
system  arising  in  the  computational  modeling  of  neuron  systems.  In  the  first  example  the  voltage 
in  a  neuron  was  modeled  using  the  highly  nonlinear  Hodgkin-Huxley  model.  The  neuron  receives 
pulses  from  three  branches.  The  output  of  interest  is  the  voltage  at  a  node  in  the  neuron.  Figure  5 
shows  an  excellent  agreement  between  the  full  order  model  with  roughly  N=1200  degrees  of  free¬ 
dom  and  the  reduced  order  mode  with  n  =  30  degrees  of  freedom.  The  application  of  the  reduced 
order  model  is  about  10  times  faster  than  the  application  of  the  full  order  model.  The  spiking 
behavior  of  the  voltage  in  the  neuron  is  captured  by  the  reduced  order  model. 
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Figure  5:  Full  and  reduced  order  modeling  of  the  voltage  in  a  neuron  modeled  by  the  nonlinear 
Hodgkin-Huxley  equations.  The  spiking  behavior  is  captured  by  the  reduced  order  model,  which 
can  be  executed  100  times  faster  than  the  full  order  model. 

Another  test  uses  the  FitzHough-Nagumo  equations  to  model  a  neuron.  Figure  6  shows  that 
the  reduced  order  model  can  capture  the  limit  cycle  behavior  of  the  original  full  order  model. 


Figure  6:  Full  and  reduced  order  model  of  the  FitzHough-Nagumo  equations.  The  reduced  order 
model  captures  the  limit  cycle  behavior  of  the  original  full  order  model. 
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